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Interference effects in the counting statistics of electron transfers through a double 

quantum dot 

Sven WelackQ Massimiliano EspositoQ Upendra Harbola, and Shaul Mukamel§ 

Department of Chemistry, University of California, Irvine, California 92697, USA. 

We investigate the effect of quantum interferences and Coulomb interaction on the counting 
statistics of electrons crossing a double quantum dot in a parallel geometry using a generating 
function technique based on a quantum master equation approach. The skewness and the average 
residence time of electrons in the dots are shown to be the quantities most sensitive to interferences 
and Coulomb coupling. The joint probabilities of consecutive electron transfer processes show 
characteristic temporal oscillations due to interference. The steady-state fluctuation theorem which 
predicts a universal connection between the number of forward and backward transfer events is 
shown to hold even in the presence of Coulomb coupling and interference. 



I. INTRODUCTION 

Fast and sensitive charge detectors and highly-stable 
current bias sources have made it possible to measure 
individual electrons crossing arrays of tunnel junctions 1 
or quantum dots 2 -! 3 -! 4 - (QDs). Directional forward and 
reverse counting through two quantum dots in a se- 
ries has been reported^. Spurred by this experi- 
mental progress, electron counting statistics (ECS) in 
nanosystems has attracted recent theoretical interest 
both in the noninteractin g 6 ' 7 ' 8 ' 9 ' 10 ' 11 ^ ^^^^ and 
Coulomb blockade regim o 17 ! 18 . It has been shown 
that strong Coulomb interactions suppress large current 
fluctuations^ 7 -. 

Double quantum dot (DQD) systems in paral- 
lel geometr y 19 ' 20 i 21 1 22 and single multi-level quantum 
dot o 23 i 24 can display interference effects due to the mul- 
tiple paths that electrons take to cross the junction 
when Coulomb interaction is present. These effects on 
the average electric current and populatio n 19 i 20 i 21 i 22 i 23 i 24 
do not require a magnetic field as in double quantum 
dot Aharonov-Bohm interferometer o 25 i 26 i 27 i 28 ' 29 i 30 . A re- 
cent study of interference effects on the electron transfer 
statistics has been shown that they can induce super- 
Poisson shot noise in Coulomb blockade regime and the 
high bias limit 18 . Using the terminology ofi^,, we denote 
the couplings between populations and coherences in the 
many-body eigenbasis as interferences. 

In this paper we extend full counting statistics to arbi- 
trary Coulomb coupling in a DQD junction by employ- 
ing a quantum master equation approach (QME). We 
find that the third moment and the average electron res- 
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idence time in the system provide useful indicators for 
interferences effects on the ECS. We also introduce joint 
elementary probabilities to electron transport which can 
be obtained from a time-series analysis of single trans- 
fer events. These probabilities reveal temporal oscilla- 
tions induced by interferences. Their amplitude is am- 
plified by the Coulomb coupling and the frequency is de- 
termined by the energy detuning between the orbitals 
of the dots. This method can reveal detunings much 
smaller than fcf,T. The employed elementary probabili- 
ties are equivalent to waiting time distributions recently 
investigated in single electron transport 3 - 1 - and electron 
transport through single molecules^ 2 -. 

ECS experiments have been proposed recently^ as good 
candidates to test the validity in the quantum domain of 
far-from-equilibrium fluctuation relations, called fluctua- 
tion theorems (FTs) , which have raised attention in clas- 
sical system o 33 ! 34 ' 35 ! 36 ' 37 . In Ref.— noninteracting elec- 
trons were considered and interferences were absent. We 
demonstrate in this paper that the FT still holds in the 
presence of interferences and Coulomb repulsion. 

Counting statistics in systems with interferences can 
be realized experimentally by connecting a DQD to two 
leads. The transitions between states can be measured 
by employing a quantum point contact 2 -! 3 -^. As a refer- 
ence point for the DQD junction, we also calculate the 
corresponding values for a QD with a single orbital where 
interferences are not possible. 

The paper is organized as follows: We present the DQD 
and QD models in sections (|I Aj) and (|IBj) . respectively. 
Equations of motion (EOM) for both models are derived 
in section (|TTJ . We then introduce the generating function 
(GF) in section ([111]). In section (JVJ) we derive elemen- 
tary probabilities for different electron detector configu- 
rations. The numerical results are presented in (|VI[) and 
we summarize and conclude in section (|VII[) . 



A. Model A: Double quantum dot 

Model A consists of a DQD connected in parallel to 
two leads. Each QD contains a single spin orbital which 
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FIG. 1: Model A: Double quantum dot injunction, each quan- 
tum dot has a single spin orbital with energy E s , connected 
to a left and right lead. The left and right coupling elements 
are equal and also T\ q = T^q, only the Fermi energies of the 
left and right lead are different, thus defining a bias voltage 
V = Epj — Ef.v U is the Coulomb interaction parameter 
between the two dots. 



is connected to the leads as shown in Fig.Q] Experi- 
mentally it would correspond to spin polarized leads in 
an infinite magnetic field where only one spin state is 
accessible by electronic excitations. This spinless DQD 
model has been used to study the effect of interferences 
on average currents withou t 20 ' 21 i 22 and wit h 25 i 26 i 27 i 30 an 
additional Aharonov-Bohm phase. 
The Hamiltonian 



H = H S + H R + H, 



SR, 



(1) 



represents the system Hs, the lead Hr, and their cou- 
pling term Hsr- We label the two dots with index 
s = 1,2. In second quantization, the local-basis Hamil- 
tonian of the system reads 



2 

E 



(2) 



U is the Coulomb coupling strength. The environment 
consists of two independent leads in thermal equilibrium. 
For each lead, the Hamiltonian is given by 



H, 



(3) 



with and *£ q create and annihilate an electron in lead 
mode \q) with energy u) q . To keep the notation simple, we 
will only refer to the left lead in further derivations. The 
right lead will be added to the hnal expressions. Since 
the leads are in thermal equilibrium, their occupation 
numbers are determined by Fermi-Dirac statistics 



= n F (LJ q -E F )5 q 



(4) 



where np(u) — l/{eP" + 1) is the Fermi function, [3 = 
1/kT, and Ep the Fermi energy. We denote the trace 




FIG. 2: Model B: A single orbital QD with energies Eif, E21 
for the spin-up and spin-down state. The coupled leads have 
Fermi energies Ef,i and Ef,t- U is the Coulomb coupling 
parameter. The left and right coupling elements are equal. 



over the lead degrees of freedom by (-)r = ^r{-Pr} 
where pr is the density operator of the lead. The cou- 
pling of the lead with the system can be written as 



sq 



(5) 



T sq is the coupling strength between lead mode q and the 
s-th QD. We assume weak system- lead coupling and no 
direct overlap (tunneling) between the wave functions of 
the left and right lead. 



B. Model B: Single quantum dot 

The second model, shown in Fig. [21 consists of a QD 
with a single orbital which can accommodate two elec- 
trons with opposite spins coupled to two leads. In anal- 
ogy to Model A, the total Hamiltonian can be written as 
H = Hs + Hr + Hsr- The system part is given by 

H s - E ^*L*i<r + ^I-r^iT^u^U- ( 6 ) 



Here, the spin is denoted by a =|, j. The Hamiltonian 
of the left lead reads 



(7) 



Electrons with different spins in lead mode \q) have the 
same mode energy uj q . The coupling term between the 
system and the left lead is given by 



H S r = £(V, ff *L*, + 



(8) 



V qa is the coupling strength between lead mode q and 
spin orbital a. Even though both Models A and B can 
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accommodate two electrons, there is a qualitative differ- 
ence between the two. In Model B the lead operators 
have spin indices so that 

= n F (uj q - £ F )<V<W- (9) 

The additional b aa i in Eq. © compared to Eq. (U|) will 
be critical, as will be shown later on. 



II. EQUATIONS OF MOTION FOR THE 
REDUCED SYSTEM DENSITY MATRIX 

Master equations have been widely used to 
simulate electron transport through quantum 
system s 38 ' 39 ' 40 ' 41 ' 42 ' 43 . Using a numerical decompo- 
sition of the lead spectral densit y 44 i 45 , a non-Markovian 
master equation for electron transport was derived 
for non-interacting 4 ^ 4 ^ and interacting^ 8 , electrons to 
second-order in H$r. Higher order coupling elements 
can be derived via path-integral calculus 3 ^,. The ap- 
proach is valid for arbitrary Coulomb coupling strength, 
temperature and bias. Master equations can be used to 
calculate the GF of the charge transfer statistics ^ 17 ' 18 . 

The total density operator is denoted p, and ps = 
tift{p} and pa = trs{p} denote the system and the leads 
components of p, respectively. We use the full Fock space 
as a basis for the system part and hereafter we will re- 
fer to ps simply as density matrix. The quantum mas- 
ter equation to second-order perturbation theory in Hsr 
reads 



Ps(t) = -iC s p s {t) 



(10) 



The auxiliary operators in Eq. are given by 

=E fMc£){t-t>)Gs{t,t')* s ,, (12) 

*S -) (*) = £ fdt'(c[-){t-t')XG s {t,t')^ s ,. (13) 

The correlation functions C^,(t) are discussed in more 
detail in Appendix JB|. Because we assume that the 
coupling of the lead to the two dots is symmetric T q \ = 
T q 2, the correlation functions (Eqs. (|B3|B4|) ) which cor- 
respond the cross coupling terms, s =/= s', can be written 
as C ( s f,\t - t') = Cif\t - t'). Setting C ss >{t) = for 
s 7^ s' would correspond to the rotating wave approxi- 
mation (RWA) 3 ^. 

Liouville-space allows a compact super operator no- 
tation. We define the following Liouville space super 
operators LA = [H S ,A], $( fl >4 = A^, = 
yKR)A = and &(L) A = ^t^. a is an arbitrary 
operator in the space of the system. We use R, L to de- 
note left and right super operators and l,r as indices for 
the left and right lead. Eq. (fTTj) finally reads 

Ps(t) = W A (t) Ps (t) (14) 

where WA(t) is 

W A (t) = -%L - H'(t) + E l + (t) + Ei_(t). (15) 

The dissipative term is separated into a diagonal contri- 
bution 



Jto s 



Eq. (|10p is derived in Appendix [X] The Liouville oper- 
ators are defined via £sr: — [Hsr, ■], £-s- — [Hg, •] and 
the propagators are Gs+R(t, t') — exp(— i(Cs + Cn)(t — 
t')) and G s (t,t') = exp(-i£ s (t-t')). We also set h=l. 
In order to propagate Eq. (fTTj)) . we derive EOM for its 
dissipative part (second term in Eq. (I10|l). This will be 
done in the following two subsections. 



which leave the number of electrons in the system un- 
changed and two off-diagonal parts 



elw = ^*i (ji) *^ i) (()+f[ L, $i (+ '' i) ((), (is) 



A. Quantum master equation for Model A 

Applying the system-lead coupling term of Model 
A to Eq. QI| one gets 4 ^ 

Ps (t) = -i[H s , Ps (t)] 

s 

+E (*i +) (*)ps(t)*i - PsWPml) 

s 

+h.c.} (11) 



which increase or decrease the number of electrons, re- 
spectively. WA(t) is a 16 x 16 matrix in Liouville space. 
The terms for the right lead n r (t), T,^(t) and £l(t) can 

be derived by replacing C^) in Eqs. (|12|13p with the cor- 
relation functions of the right lead and have to be added 
to Eq. 1(15]) ■ 

In Appendix (|B]) we discuss the spectral decomposi- 
tion of the correlation functions into the form (t) = 

Y^k=™ a i ± ^ e7 * Bi-directional counting requires rel- 
atively small bias voltages in order to have significant 
backwards transfer rates. The Fermi function is ex- 
panded into Matsubara-frequencies, i"f^ for k = m . + 
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1, m', and small bias voltages can be used in our cal- 
culations. A numerically efficient way for calculating 
Eqs. (|12|13p is by propagating the EOMi& 

+E«i ±} ^- (19) 

s' 

Summation over the spectral decomposition given by 
Eq. (|B10[) results in an explicit expressions for the many- 
body auxiliary '(t) — YHk=\ % fc (*) operators. 
Eqs. (TTJ] - [T9"l) propagated simultaneously in time form 
a closed set for calculating the reduced density opera- 
tor ps(t). An adiabatic switching of the coupling is not 
required in this method. We are interested in electron 
counting statistics at steady state which we reach by nu- 
merically propagating pg — ps(t — > oo), II = H(t — > oo), 
S + = — > oo) and E_ = S_(f — > oo). 

The KWA as applied in Ref.?- decouples the coherence 
and population part of the system density matrix in Li- 
ouville space eliminating all interferences. For Model A, 
the RWA applied to Eq. (fT4"]) would reduce it to a Pauli 
rate equation in the eigenbasis of the system. This formal 
cancellation can be realized as is illustrated by Model B, 
a single QD. Model B also allows a maximum number 
of two electrons, but interferences are absent due to the 
different spin states of the electrons. In the following 
subsection, we present the Pauli rate equation for Model 
B without invoking the RWA. 

B. Pauli rate equation for Model B 

Inserting the system-lead coupling term ([6]) of Model 
B into Eq. ([5]), we can write 

Ps (t) = -i[H s , Ps (t)} 

a 

+h.c.} (20) 
The auxiliary operators in Eq. (|20|) are given by 

¥+J (t) = f dt'c£> (t t')G(t, t')*i„, , (21) 

J to 

*£>(i) = fdt'(c^{t-t')XG{t,t')^ la ,. (22) 

J to 

The correlation functions ci^\t — t') are given in 
Eqs. (|B1IB2[) . The master equation can be recast in the 
compact Liouville space form 

PS® = W B (t) Ps (t) (23) 



with 

W B (t) = -iC - U l (t) + E l + (t) + E l _(t). (24) 
The diagonal part is given by 

m = E^ L3 ^' i) (*)+*Il R) M;' R) (*)+/ l -c. (25) 

and off-diagonal parts are 

Zi(t)=^\?H? R \t) + *[?*ll~' L) (t), (26) 

a 

ELM = E *l R) *£' L) (t) + *\M?' R) (t)- (27) 

Because Cj^ = Cyt' — (see Appendix ifH])) the spin 
quantum number causes a separation of coherence and 
population part in the density operator. The population 
part of Eq. (|2"3")l satisfies the Pauli rate equation 

P s (t) = W P P s (t). (28) 

P denotes the population of the states and Wp is the 
Pauli rate matrix. 

III. THE GENERATING FUNCTION 

The transfer probability of k = (fcijfea) electrons in 
the time interval t — to through the left lead is denoted 
P(k:t). We define the generating function G(\;t) by 

G(i\;t) = EX fc Mfc;i) (29) 

k 

The probability distribution is obtained using P(k;t) = 

&J?*G(iA;t)e-** 

We will calculate the GF by tracing the generating 
operatop£ii£ (GO) <?(A; t) over the system degrees of free- 
dom 

G(\;t)=tx s {g(\;t)}. (30) 

Based on the master equation (|2"3"|) . g(X;t) satisfies the 
EOM 

g(\;t) = W(X)g(X;t). (31) 

Counting starts after the system has reached steady 
state. Thus the initial condition of the GO is given by 
the steady state of the system density matrix g(X;t = 
0) = pf. The propagator W{X) of the GO is 

W{X) = -i£-n l -n r +e Xl Y, l + + e X2 Y, l _+Y, r + + Y, r _. (32) 

We shall consider the statistics of charge transfers be- 
tween the left lead and the system. A := (Ax, A2) controls 
the specific statistics obtained by propagating Eq. (|31[) . 
We investigate 3 cases. Setting A + := (Ax = A, A2 = 0) 
gives the counting statistics of the incoming electrons, 
A_ := (Ax = 0, A2 = A) the outgoing, and A„ := (Ax = 
A, A2 = — A) the net-process. The corresponding proba- 
bilities are denoted by P^(k; t) with rj = +, i, n, respec- 
tively. 
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IV. CUMULANTS OF THE TRANSFER 
DISTRIBUTIONS 

We shall calculate the first C\(t), the second Cz(t) and 
the third C^t) cumulant of P^ v )(k;t) with respect to 
k. Ci(t) = k/t = Y,k kP ( k '^ t )/ t is related to the av- 
erage current I(t) = eCi(t). The second cumulant de- 
fined by C^t) = (k 2 — k 2 )/t is commonly represented 
by the Fano factor F(t) = q'Ya ■ The third cumulant 



Cs{t) = (k — k) 3 /t measures the skewness of the prob- 
ability distribution P(k; t) with respect to k. For small 
bias, the current is small and the electrons transferring 
into the system are uncorrelated. This leads to a Poisson 
counting statistics where the Fano factor equals F = 1. 
If F < 1, the process is Sub-Poissonian, while F > 1 
indicates Super-Poissonian statistics. 

The time-dependent cumulants can be calculated from 



Cf\t) = d{K{\ n ;t)\ x=Q 
where K(X; t) is obtained from the GF 



(33) 



(34) 



We make use of the fact, that K(X; t) has a well defined 
long time limit and calculate the asymptotic values^ for 
the steady state current 



C[ v) = — lim K(X n ;t) 

OA t— >0O 



A=0 



and the zero-frequency power spectrum 



d 2 



(35) 



(36) 



A=0 



The asymptotic value of the skewness is defined in the 
same fashion 



a 



(</> 



Q3 

— l}mK(X v ;t) 



(37) 



A=0 



The long time limit ofK(X v ; t) can be calculated from the 
dominant eigenvalue^ ei(A^) of the propagator W(X V ) 
given by Eq. Then, K(X V ) = lim^oo K(X n ; t) = 

e^Xr,). Note that e^A,, = 0) = 0. 



V. PROBABILITIES OF ELEMENTARY 
EVENTS 

We introduce elementary probabilities to characterize 
consecutive electron transfer events m\ 1 m2 at times ti,t?,. 
m characterizes the side of the process (I, r) and if the 
electron is transferred in (+) or out (-) of the QD. The 
elementary probability is given 



P(t 2 ,tt 



l\pm 2 c v m i Q \ 



(38) 



Eq. (|38|) is the joint probability of detecting specified elec- 
tron transfers at times t\ and ti with no transfers oc- 
curring in the time intervals ti,to an d t^t\. We de- 
note the trace over the system degrees of freedom by 
(•} = trs{-ps(io)}. Sti.tj is the propagator of the system 
in absence of transfer events at the leads within the time 



interval ti, tf 



exp ((-iC IT) {U - tj)) 



(39) 



We shall consider three cases, i) An electron is de- 
tected when it enters the junction through the left lead 
at time t\ = and leaves the junction through the right 
lead at time t. The electron transfer operators are H\ 
and £2 respectively and we denote this transfer pathway 
by I — > r. The joint elementary probability is then given 

by 



fUr(Mo) = (ElSt.toSL). 



(40) 



ii) The reverse process, r — > I, which we write as 

^z(Mo) = <£-Smo£;>- (41) 

Hi) Transfer from the left lead into the system and back 
into the left lead denoted by Z — ^ Z can be written as 



(42) 



Quantities (|40j|42p can be measured as follows. One 
has to detect single directionally resolved electron trans- 
fers between the leads and the system and record a suf- 
ficiently long time-series of transfer events. Then a his- 
togram of the number of occurrences of a specific con- 
secutive transfer event, such as I — > r, as function of 
increasing time intervals t — to can be generated. The 
histogram has to be normalized by the total number of 
events S mi — > S m2 in the time series. 

The total conditional probability P^ ll ^ ni2 to have an 
S™ 2 event if the last event was a £ mi event irrespective 
of the time-interval between 7712 and mi can be calculated 
by integrating the time-dependent conditional probabil- 
ity P^ m M^) = (£™ 2 5 Ti t £ mi )/ (£ mi ): 



P' 



dr(E m2 5 r , to I] mi )/(S mi ). 



(43) 



Other interesting quantities can be calculated from 
Eqs. (|40j(4"2"l) . For example, the residence time of elec- 
trons in the system subject to a specific transfer process 
rni — » 7712 is given by 



/~dr<E"*<>S r , t0 E"*i) ' 



(44) 



We propose another interesting setup by counting elec- 
trons only at the left lead regardless of electron transfers 
occurring at the right lead. Thus the propagator ([39f is 
modified to 



S Ulti = exp ((-»£ -tf - IT + S; + SI) (ti - tj)) 



(45) 
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The conditional probability of an electron entering the 
system at time t Q through the left lead, and the next 
electron entering at time t from the left lead can be writ- 
ten as 

P^(t\t ) = + S t<to E l + )/(E l + ). (46) 

We denote this transfer series by I, I. In the time in- 
terval t-t electron transfers take place only at the right 
lead. For uncorrelated electron transfer, the propability 
distribution ([46)) is Poissonian and can be written as^ 

P% i ''(t\t ) = e- c «(I* + ). (47) 



8 3 = E^+U 



81=^=0.999 



FIG. 3: Level scheme in many-body eigenbasis for Model A. 
For Model B, E\ and E2 have to be replaced by and 



VI. NUMERICAL SIMULATIONS 

The master equation for Model A ([H]) and the EOM 
of the auxiliary creation operator (I19|) form a system of 
equations of motion and we propagate them simultane- 
ously into steady state using the Runge-Kutta method. 
At steady state, the memory of the non-Markovian mas- 
ter equation vanishes and Eq. (fT4"|) corresponds to its 
Markovian counterpart and the electron transfer oper- 
ators Eqs. (|16H18J) become time- independent. Using the 
obtained transfer operators and steady state density ma- 
trix, we then propagate the EOM of the generating op- 
erator (|3Tj) to a finite binning time t for each step of the 
discrctized counting field A to derive the time-dependent 
cumulants. Alternatively, we use an eigenvalue decompo- 
sition of the steady state propagator of Eq. (f52"| in order 
to calculate the asymptotic cumulants. The same ap- 
proach is used for Model B based on Eq. (|2"5)) and its 
dependencies. 

We shall expand all system operators in the many- 
body eigenbasis of the system. For Model A, we use the 
transformation 

mm' 

The many body basis is spanned by four states |0) = |00), 
1) = |01), \2) = |10) and |3) = |11). The coefficients 

(t s) 

of the transformation a^'^, between orbital basis and 
many-body basis can be derived from the Fermion anti- 
commutator relations. Here the non-zero coefficients for 
the creation operators are o^q = 1, 0$ — 1- aiyj = 1 
and a§ = -1. Thus we get *| = |10)(00| + |11)(01| and 
= |01)(00| - 1 11) <10|. The annihilation operators can 

be derived by replacing a^ m , with or^, m . For Model B, 
one has to replace s — 1 with a —\ and s = 2 with =[. 
We can write @ and ^ as 

3 

H s =^e m \m){m\. (49) 
m=0 

e m is the energy of state \m). The energy scheme we 
used for Model A and B is shown in Fig. ©■ All ener- 
gies are scaled with respect to the equilibrium chemical 



potential of the leads E F = 1. The fixed orbital energies 
for Model A are e x = E 1 = 0.999, e 2 = E % = 1.001 and 
for Model B e x = E^ = 0.999 and ei = E n = 1.001. 
Thus, the double occupancy state has an energy of £3 — 
Ei + E 2 + U = 2.0 + U. The ground state is set to zero 
60 = 0. The bias voltage is applied symmetrically to the 
system using E F ,i = E F + V/2, E F . r = E F — V/2. To 
study temperature fluctuations of electron transfer as a 
function of U, we set T = 0.002 to be in a temperature 
range of (3V ~ U . In Eq. (|B6|) . we restrict the spec- 
trum to a single Lorentzian centered at fix = E F . We 
also chose relatively large bandwith parameter Ti = 1 
(wide-band limit). We set p 1 = 2 • 10~ 4 in Eq. (|B6|) 
(pi ~ ^2 q l^?| 2 = J2 q \T q \ 2 ) as a reasonable small value 
considering the weak coupling required in order to guar- 
antee physical results within second-order perturbation 
theory. 



A. The bidirectional transfer probability 

In Fig. 0] we show the probability P^ n '{k\ t) of the net- 
number k = k\ — ki of electron transfers between the 
left lead and the system at steady state. The binning 
time is fixed to t — 200. The left panel depicts Model 
A, the right panel Model B. The distributions for differ- 
ent Coulomb coupling strengths U = , U = 2 • 10~ 3 
and U — 8 ■ 10~ 3 are fitted with a Gaussian P(k) = 

-0Lexp(—^(k—x c ) 2 ). The fit parameters are given in the 

graph. A small bias of V — 4- 10~ 4 allows electron trans- 
fers against the bias. As discussed for non-interacting 6 
and strongly coupled^ electrons, the Gaussian provides 
a good approximation for P^ n '[k\t). This is mainly due 
to small non-equilibrium contrains. For high bias, the 
deviations will be more significant. Interference effects 
are enhanced by Coulomb interaction. This can be seen 
by comparing the results for Models A and B in Fig.0] 
where interferences are present and absent, respectively. 
For U = 0.0, the probability distributions for both mod- 
els are the same. With increasing U, the fitted curve is 
shifted and broadended in the presence of interferences. 
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FIG. 4: Electron transfer probability P n (k,t) for different 
Coulomb coupling strengths U as function of the net-number 
of transferred electrons k. Left panel depicts Model A, right 
frame Model B. The time is set to t — 200. A small bias 
is applied V = 4 • 10 -4 . w —variance and x c —position of 
peak are the fit parameters of the Gaussian. In Model B, the 
curves for U = 2 • 10" 3 and U = 8 ■ 10~ 3 overlap. 
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FIG. 5: Cumulant generating function K(t,X n ) of the net- 
process as a function of A for different values of Coulomb 
repulsion U. Left panel shows results obtained for infinite 
binning time limt^oo K(t, A n ), right panel for a finite binning 
time of t — 10. The x-axis is rescaled by [5V. 



Cumulants of the transfer probability 
distribution 



B. The fluctuation theorem 



Fig.E] shows K(t, A n ), given in Eq. (|34j) . for the net 
process as a function of A for different Coulomb cou- 
pling strengths. The x-axis was rescaled by the fixed 
bias of V — 2 • 1CP 3 . To explore the role of in- 
terferences we compare Models A and B. The left 
panel shows the asymptotic cumulant generating func- 
tion lim^oo K{t 1 A„) which is computed from the eigen- 
values of W(X n ), the right panel shows K(t,X n ) for fi- 
nite binning time t = 10 calculated by propagation the 
EOM (|3"Tj) . We observe that finite Coulomb coupling, 
U = 2 ■ 10 ~ 3 , significantly changes the cumulant GF 
K(t,X n ). This indicates that the influence of Coulomb 
blockade can be measured in the first cumulant, the 
current. The effects due to the interference terms are 
smaller and are observable for U ^ 0. In this regime, 
a discrepancy between Model A and B can be observed 
for X/PV w 0.5 which indicates that it is significant in 
the higher cumulants. We find that at infinite binning 
time the symmetry K(t, \)=K(t, A — j3V) holds in both 
models as shown in the left panel of Fig. [51 Immedi- 
ately follows from this symmetry the fluctuation theorem, 
pr% = exp(/3Vk), for t — > oo. This relation thus holds 
in systems with interferences and Coulomb interaction. 



In the following, we discuss the effects of quantum in- 
terference and Coulomb coupling on first cumulant Ci, 
the Fano factor F = C2/C1 and the normalized skewness 
C3/C1 for a infinite binning time. 

The upper panel of Fig.[S] depicts the first cumulant 
through the left lead for three different processes: A + 
(ingoing), A (outgoing) and X n (net-process). The bias 
is V = 1 ■ 10~ 3 . U introduces an energetic penalty for 
double occupancy. This explains the current drop around 
U = 2-10~ 3 by the fact that the energy of the double 
occupancy state increases with U. This reduces the oc- 
cupancy of the state and its contribution to the total 
average current. The drop is smooth due to the finite 
temperature. Interesting effects can be observed at the 
intermediate coupling regime. So is the drop of the net- 
current accompanied with an increase of the reverse cur- 
rent. The Fano factor and the normalized skewness of the 
transfer in direction of the bias as well as the net-transfer 
show an increase following the decrease in current due to 
increasing U. The corresponding values for the reverse 
transfer show a weak response. 

Comparison of Model A with Model B shown in Fig.[7| 
demonstrates the influence of interferences on the net- 
transfer statistics (A„). The skewness is of particular 
interest since it was found to be the most sensitive of the 
three in systems with interference^. The Fano factors of 
Model A and B diverge in Fig.[7| with increasing U but 
their absolute difference remains small. The skewness 
reveals the presence of interferences more clearly than 
the Fano factor. The average current shows a negligible 
difference. The observability of interferences requires a 
strong Coulomb coupling since the differences between 
Model A and B in the third cumulants are significant 
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Model A 




U (10" 2 ) 

FIG. 6: First cumulant Ci (top panel), Fano factor C2/C1 
(middle panel) and the normalized skewness C3/C1 (bottom 
panel) as functions of Coulomb coupling U for infinite binning 
time. A+ refers to the incoming, A_ the outgoing, and A n the 
net-transfers of electrons between the left lead and the system. 
The bias voltage is V = 1 ■ 10~ 3 . 

when the double occupancy state is inaccessible by elec- 
tronic excitations. 



D. Joint elementary probabilities 

Using expressions (|40H42I) . we calculated the joint 
probabilities of directionally resolved consecutive elec- 
tron transfers. Fig. [5] compares the joint probabilities 
P/_ r (i, t ), Pi^iit, t ), P r ^i(t, to) for the two models and 
a small bias of V = 2 • 10~ 3 . In Model B, we observe an 
almost exponential decay of the probability of an out- 
going electron transfer event at the left/left /right site 
following an incoming electron transfer event at to at the 
right/left/left site. A weak non-exponential slope indi- 
cates a weak correlation between the transfer processes. 
The probabilities of the reverse processes r —> I decay 
faster than in the direction of the bias / — > r. Conse- 
quently, the I — > I process has an intermediate decay 
rate. Electron-electron coupling leads to a slower decay 
in all three processes, indicating that the probability of 
an electron to reside longer within the junction increases 
with U . The probability has its maximum at to since the 
orbitals are in direct contact with both leads. 

In Model A, the exponential decay of the probability 



Model A 




U (10" ) 

FIG. 7: Comparison of first cumulant, Fano factor and 
normalized skewness for the net-transfer statistics (A n ) be- 
tween Model A and B for infinite binning time. The bias is 
V = 1 • 10" 3 . 



is superimposed by an oscillation due to orbital inter- 
ference. The amplitude of the oscillation is increased 
by Coulomb interaction. We Fourier transformed the 
P/_> r (t, to) probability of Model A in order to analyze the 
dependency of the frequency on the system parameters. 
The magnitude of F(uj) = -j= J °° dt e - luJt Pi^ r {t, t = 0) 
is shown in Fig. [9] From the peaks we can conclude that 
the frequency LOi^ r of the oscillation in Pi—, r (t,to = 0) 
is determined by the detuning of the orbitals of the two 
quantum dots u>i^ r = E\ — E% = e\ — £2. We could 
not find a dependence of wi^ r on U . The energy gap 
between the double occupancy state and the single occu- 
pancy states £3 — ei (£3 — £2) is three orders of magnitude 
larger than the gap between the single occupancy states. 
We also found that the amplitude of the oscillation de- 
creases with increasing detuning. Thus oscillations due 
to the £3 — £1 (£3 — £2) energy gap which would depend on 
U have a high freqency and a small amplitude invisible to 
our numerical method. The time-series analysis of con- 
secutive electron transfer events seems to be highly sen- 
sitive to very small energy differences when interferences 
are present. Note that this method works for energy dif- 
ferences smaller than fc^T, here £1 — £2 ~ k^T, which 
cannot be resolved by average current measurements by 
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t t 

FIG. 8: Elementary probabilities Pi^ r (t,to)/ Pi^i(t,to)/ 
P r -,i(t,to) for an electron entering the junction through 
the left/left/right lead at time to and leaving through the 
right/left/left lead at time t, respectively. The left side 
shows Model B, the right side Model A. The bias voltage 
is V = 2 ■ 1CT 3 . 



scanning the voltage of the leads. 

The upper panel of Fig.[]I)]depicts the total conditional 
probability P£ ll ^ m2 as a function of external bias voltage 
for Model A (bottom panel Model B). The two transfer 
processes under consideration are I — > r, r — ► I, and we 
use U = 0.0 and U — 2 • 10~ 3 for comparison. Since the 
probability in Eq. (|43f is conditional on the first trans- 
fer, the r — ■> I and I — > I processes are equally likely. In 
the absence of an external bias, P[ C _, r and -P°->i are ecma l 
and the electron transfers are driven by the thermal fluc- 
tuations of the leads at finite temperature. A splitting 
occurs when the voltage is turned on. The deviations 
caused by the Coulomb coupling vanish at high bias in 
both models. 

Fig-El uses the same parameters as Fig.[TUl but the 
observables are shown as functions of the Coulomb cou- 
pling U at a fixed bias of V = 2 • 10~ 3 . The upper 
panel depicts the total conditional probabilities, the bot- 
tom panel the average residence times for Model A and 
B. Increasing U causes a dip of the / — > r probabilities 
around U = 2- 10~ 3 . There are no qualitative differences 



2.5x1 () 




co (10 3 ) 

FIG. 9: Magnitude |f of the Fourier transformation of 
Pi-> r (t, to) for Model A. AE — E1—E2 is the energy difference 
between orbits of the quantum dots. Coulomb coupling is set 
to U = 4 ■ 10" 3 . The bias voltage is V = 2 • 10~ 3 . 
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FIG. 10: Total conditional probabilities (Eq|43} Pf_ tr and 
Pr-,i as a function of bias voltage for different Coulomb cou- 
pling strength U. Upper panel Model A, bottom panel Model 
B. 

between the probabilities of the two models. U strongly 
affects the mean residence time of the electrons which 
is also sensitive to the presents of orbital interferences. 
Thus the residence time between the electron entering 
and leaving the junction provides a useful measure for 
the presence of Coulomb coupling and interferences. 

Fig-EH depicts the conditional joint probability 
Pfl(t\to), Eq. ([46]) . for an electron to enter the junction 
through the left lead at time to and next electron to enter 
at time t also through the left lead. The detector is only 
applied to the left lead, transfers at the right lead are per- 
mitted at all times. Models A and B are compared for 
different U. The probabilities of Poissonian two-electron 
transfers, Eq. (|47p. are shown for comparison. Similar to 
the configuration with electron detectors at both leads, 
Fig. [5J we observe longer tails in the probability distribu- 
tion in time for larger U. Consecutive electron transfers 
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FIG. 11: Upper panel: Time-integrated conditional joint 
probabilities (Eq |43|l . P{L, T and P^-fi versus Coulomb cou- 
pling strength U. Bottom panel: The mean residence time 
t res of an electron in the system as a function of U for the dif- 
ferent processes and Model A and B. The bias is V = 1 • 10 -3 



VII. CONCLUSIONS 

We have calculated the counting statistics in a model 
junction for finite Coulomb coupling strength, bias and 
finite temperatures. The numerical results reveal several 
significant measurable effects of the quantum interference 
in a DQD on the electron transfer statistics. 

The skewness provides a sensitive observable. We also 
show that a measurement of the average residence time 
of electrons is affected by quantum mechanical interfer- 
ence as well as the Coulomb coupling between two paral- 
lel quantum dots. The observed oscillations in the joint 
elementary probabilities can be recovered by time-series 
analysis. Their amplitude and frequency are directly re- 
lated to the Coulomb coupling and the level detuning of 
the DQD, respectively. 

Several extensions of the model could be of interest. 
Decoherence effects can be included by coupling a dis- 
sipative phonon bath to the sites. Including higher or- 
der coupling elements beyond second-order perturbation 
theory in the GF could reveal additional insights into the 
dynamics of the system- lead contact. 



o.oi 



\. 1 1 1 1 


1 






Poisson, U=0.0 






Poisson, U=0.01 






Model B, U=0.0 






Model B, U=0.01 




',1,1 ~~ -i 1 ^- 1 -T 











Poisson, U=0.0 






Poisson, U=0.01 






Model A, U=0.0 






Model A, U=0.01 









t 

FIG. 12: Conditional probabilities Pfi(t\to) for an electron 
entering the junction through the left lead at time to and the 
next electron enters at t. Only transfers through the left lead 
are detected. Electrons are allowed to leave and enter through 
the right lead at all times. Model A and B are compared for 
different U. For comparison, we also plot the corresponding 
Poissonian processes P ; p ° lss (f|to). The bias voltage is V = 
2 ■ 10" 3 . 



are strongly correlated at short times leading to small 
probabilities for the next electron to enter the junction 
after the first one has entered. At intermediate times, the 
probability is larger than the probability of the Poisso- 
nian process. For long times, the transfer becomes weakly 
correlated and the probabilities for Model A and B are 
close to a Poissonian distribution. The oscillations are 
present due to orbital interference in Model A only. 
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APPENDIX A: DERIVATION OF THE 
QUANTUM MASTER EQUATION 

In this section, we present the derivation of the 
Nakajima-Zwanzig operator identit y 53 ! 54 ' 55 and utilize 
it to couple an lead to the relevant system, here the 
quantum dots. We define a projection operator P, with 
P 2 = P, which acts on an arbitrary operator A defined 
in the Hilbert space of the full system 

PA = Btv R {A}. (Al) 

B is an operator defined in the lead part only. Applying 
the projection operator to the density operator of the full 
system leads to 

Pp(t) = p R tr R {p(t)} = p R ® p s (t). (A2) 

Thus, the EOM for the full system can be written as 

p(t) = -i£p(t) = -iCPp(t) - iCQp{t) (A3) 

which can be interpreted as the evolution of the projected 
part Pp{t) plus the evolution of its orthogonal comple- 
ment Qp = (1—P)p. Applying the projection operator to 
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the Liouville equation and to its orthogonal complement 
leads to 

Pp(t) = -iPCPp{t) - iPCQp{t) (A4) 

and 

Qp(t) = -iQCP P {t) - iQ£Qp(t), (A5) 

respectively. Integrating the differential equation of the 
orthogonal complement part (|A5|) and applying it to Eq. 
(|A4j) results in the Nakajima-Zwanzig operator identity 

Pp = -iPCPpit) (A6) 
-iP£Te- lf to d ^- p ) £ (l _ P) p ( to ) 

+PC / dt'fe- t ^ dT ^- p ^(l-P)CPp(t') 

J tg 

which is valid for arbitrary time-dependent Hamiltoni- 
ans. T is the positive time-ordering operator. One can 
further simplify by tracing over the lead degrees of free- 
dom and employing the property of the projection oper- 
ator tr n{Pp(t)} — ps(t). This gives 

Ps{t) = -i£sps{t) - itTR{C S RPii} 

+ [ dt'K(t,t') Ps (t') + In(t) (A7) 
J to 

where the Kernel K(t,t') reads 

K{tJ) = -trRiCsiife-'tiW-W 

x(l-P)(£ R + £ SR ) PR } (A8) 

and the initial value term is given by 
In(<) = -itv R {C s fe~ 1 ^ d ^- p > £ (l - P)p(t )}. (A9) 

In order to derive a practical method for solving Eq. (|A7|) , 
we utilize second order perturbation theory by assuming 
that (1 — P)Csr ~ or Csr ~ P£sr, what implies for 
the Kernel K(t, t') that 

e -i{l-P)Ct _ e -i(l- P){Cs+C R +C S R)t K e -i(l-P)(C s +C R )t 

(A10) 

We can further simplify Eq. (|A7|) by making use of the 
following exact relations: C R p R = 0, PCsrPr = 0, 
e -i{\-p)Ct = p + (i _ p) e ~ ia , tT R {C SR PC} = for 
an arbitrary operator C. Finally, neglecting the initial 
value term (|A9[) in the evolution of the reduced density 
matrix results in the Liouville equation with the system- 
lead coupling term in a time-nonlocal (TNL) regime 

Ps (t) = -iC s ps{t) (All) 
~tT R {C S R ( dt'Gs + R(t,t')CsRp(t')} 

Jt 

We neglect the influence of dissipation during the 
propagation of the density operator by applying the 
substitution^ 

P (t') = Gl(t,t')p(t) (A12) 



to the TNL Kernel. This leads to the time-local (TL) 
description given by Eq. (|10[) . For the given system, both 
TNL and TL description produce similar results in the 
weak coupling regime^,. 

APPENDIX B: NUMERICAL DECOMPOSITION 
OF THE SPECTRAL DENSITY 

In Model B, the trace over the lead degrees of freedom 
in the dissipation term of Eq. (|10p can be recast in terms 
of correlation functions of the form 

^(^Ev^^^'V^ (bi) 

cizhQ^^KA^**^****^*)* ( B2 ) 

The properties of the trace lead to C\f = c{f = 0. 
While not an approximation, the procedure mimics a ro- 
tating wave approximation in Eq. (|23p . The correlation 
functions in Model A used in Eqs. (| 1 21 13[) are given by 

C%\t) =J2T q sT; s ,(^ q e- lHRt ^ HRt pR)R (B3) 

9 

C { J(t) =J2WAe- iHRt *te iHRt %P R U (B4) 
i 

The coupling to the DQD is assumed to be symmetric 
T q i = T q 2- The cross coupling correlation functions are 

given by C^\t-t ') = cif \t-t'). The following deriva- 
tions refer to Eqs. (|B1IB2|> . The same procedures have to 
be applied to Eqs. (|B3IB4[) as well but we will not present 
them in detail. All the external properties of the lead are 
described by a single quantity, namely the spectral den- 
sity Jr(uj), which can be generated by a superposition of 
weighted 5-functions 

J R (u)=^-K\V q \ 2 5{u-u q ). (B5) 

9 

We apply a numerical decomposition of the spectral den- 
sity to derive equations of motion 

m p 2 

W =^p fc ^ (B6) 

fc=l v K > ' k 

With the complex roots of the Fermi function and of 
function (|B6|) . the theorem of residues applied to Eqs. 
ABU lB2l) results ini£ 

m 

£7<+>(i) - ^ Pfc r fc (n i ,(-^T+^)e-^ t ) 
fc=i 

-^T, J ^ e ~ w ' kt ( B? ) 
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m 



k=l 

2i 

'J 



iu k t 



(B8) 



for a =1 and a =|. We use the abbreviations = 
Qk + iTjt, = ^ — i^k and the Matsubara frequencies 
Uk given by V], — i 2wk + w -f In general, one has to 
take an infinite number of Matsubara frequencies into ac- 
count, but it was demonstrated that the summation can 
be truncatedii^^. From Eqs. (|B7IB8|) . we can write 



the correlation functions as 



m+m 



C£Ht)= E ^e^- (B9) 
fe=i 

The same decomposition can be derived for C$(t): 

m+m 

C£\t)= E a^erF*. (BIO) 



fc=i 
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